Spatially inhomogeneous phase in the two-dimensional repulsive Hubbard model 



OO 

o 
o 

(N 
o 

o 

o 



I 



I 

o 



> 
m 

OO 

in 
o 

OO 

o 



Chia-Chcn Chang and Shiwei Zhang 
Department of Physics, College of William and Mary, Williamsburg, VA 23187 

Using recent advances in auxiliary-field quantum Monte Carlo techniques and the phaseless ap- 
proximation to control the sign/phase problem, we determine the equation of state in the ground 
state of the two-dimensional repulsive single-band Hubbard model at intermediate interactions. 
Shell effects are eliminated and finite-size effects are greatly reduced by boundary condition inte- 
gration. Spin-spin correlation functions and structure factors are also calculated. In lattice sizes 
up to 16 X 16, the results show signal for phase-separation. Upon doping, the system separates 
into one phase of density n = 1 (hole-free) and the other at density ric (~ 0.9). The long-range 
antiferromagnetic order is coupled to this process, and is lost below ric- 

PACS numbers: 71.10.Fd, 02.70.Ss 



I. INTRODUCTION 

The Hubbard model^ provides a minimal framework 
for describing electron interactions in a crystal lattice, 
and has played a central role in condensed matter and 
quantum many-body physics. Especially since the dis- 
covery of high- Tc superconductors, the two-dimensional 
(2-D) Hubbard model, believed to contain the essential 
physics of the CuO planed, has been intensely studied. 
The combination of theoretical and numerical techniques 
has made important progress^, but some basic questions 
have remained. 

One of the questions is whether there is phase separa- 
tion (PS) in the ground state of the Hubbard model. The 
question is important in its own right, as a key element 
in our understanding of the phase diagram of this funda- 
mental model. Recent experimental indication of spatial 
inhomogeneities in cuprates^ has further increased its po- 
tential relevance and interest. In the past two decades a 
large body of numerical work has been devoted to resolv- 
ing this issue^'i'^'^'^iii'^ii^'^, but the results have been 
conflicting. The differing answers underscore the chal- 
lenges; the requirement of high accuracy, as well as the 
difficulty in extrapolating to the thermodynamic limit 
because of extreme sensitivity of the signal to both finite- 
size and shell effects. 

In this paper, we apply recent advances in auxiliary- 
field quantum Monte Carlo (QMC) technique a^^i^^ to 
study the ground state of the repulsive 2-D Hubbard 
model. Our goal was to shed light on the question of 
PS. A second motivation comes from ultra-cold atoms, 
where rapid experimental progress promises a new av- 
enue — optical-lattice emulatorsii — for direct "simu- 
lations" to investigate properties of Hubbard-like mod- 
els. Detailed, accurate numerical data would allow quan- 
titative benchmark and comparisons in future optical- 
lattice experiments. In our approach, the ability to con- 
trol the sign/phase problem with a good approximation, 
combined with a boundary condition integration tech- 
nique, drastically reduces the finite size and shell effects. 
This allows us to reach much higher accuracy than pre- 
viously possible in the model. The measured equation 



of state and spin-spin correlations, in lattice sizes up to 
16 X 16, show clear signals for PS at intermediate inter- 
action strengths. The nature of this spatially inhomoge- 
neous state is examined. 

The Hamiltonian for the one-band Hubbard model is: 
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-t 



h.c. 



(1) 



where cj^ (cj,o-) creates (annihilates) an electron with 
spin cr (cr =ti J.) at lattice site j, and 5 connects two 
nearest-neighbor sites. The square lattice has size ~ 
Lx L, with spin-(T electrons. The model has only two 
parameters, the strength of the interaction U/t (we will 
set t = 1) and the electron density n = (7V| + Ni)/N. 

PS occurs when the stability condition d^e{n)/dn^ > 
is violated, where e(n) is the ground-state energy (per 
site) at density n. The critical value of n can be identifled 
by Maxwell construction. Emery et al^ showed that in 
the Hubbard (or t-J) model one could study 
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e(l-/i)-e(l) 
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where h is the hole density: h = 1 — n. If PS exists, there 
is a minimum in eh{h) at he (or in the thermodynamic 
limit, a constant eh{h) for h < hc)r^- 



II. METHOD 

A. Twist-Averaged Boundary Condition (TABC) 

The signal for PS from Eq. ^ requires the slope of the 
equation of state, i.e., accurate numerical determination 
of small energy differences in the region where h is small. 
For a finite lattice, the shape of the Fermi surface varies 
considerably with n, which causes large variations in the 
energy. For example, with the usual periodic boundary 
condition (PBC), the smallest h accessible by a closed- 
shell system is ~ 0.15 in a 16 x 16 latticei^; even at 40 x 40 
the finite-size effect is still sizable, especially in the region 
relavant for PS (see inset in Fig. [31). To reduce shell 
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FIG. 1: (Color online) Upper panel: Ground state energy 
per site e(n), versus density, of the 3x3 Hubbard lattice at 
[7 = 4 (blue) and 8 (red) calculated by ED (empty symbols) 
and our QMC method (filled symbols). At each density, the 
result is the average from 1000 random © values and the 
statistical error is estimated from their distribution. Bottom 
panel: Relative error (see text) of QMC ground state energy 
compared to the exact result (percentage). 



and finite-size effects, we use twist-averaged boundary 
condition (TABC) ^^'^°i^^ , under which the wave function 
^'(ri, r2, . . .) gains a phase when electrons hop around 
lattice boundaries: 



(3) 



where L is the unit vector along L, and the twist angle 
= {9x,9y) is a parameter. With a generic 0, there 
will be no degeneracy in the one-electron energy levels. 
We average the results over many random twist angles^! 
in each system for convergence. As shown in Figs. [2] 
and [31 TABC essentially eliminates any shell effect. The 
disadvantage is that it turns the QMC sign problcnJi 
into a phase problemi^. 



B. Constrained Path Monte Carlo under TABC 

To treat this problem, we extend the constrained path 
Monte Carlo (CPMC) methodic to a Hamiltonian under 
TABC. For each given system (specified by N, n, U, and 
0), the method obtains a Monte Carlo (MC) representa- 
tion of the many-body ground state I^^g), by importance- 
sampled branching random walks (RWs) ^^i^^ in the space 
of Slater determinant wave functions. The usual sign 
problem under PBC is caused by the symmetryi^i^ be- 
tween a Slater determinant \(f)) and a degenerate part- 
ner — \(f>) (exchanging two orbitals). To specify \'^g), we 
need either, but not both. It can be show n^^i^^ that con- 
straining the RWs to {'^G\(j)) > is an exact boundary 
condition that eliminates the sign problem. In the con- 
strained path approximation, a trial wave function |5't) 
is used in place of |^'g)- 



Under TABC, the Slater determinants become com- 
plex, and we need to break the phase symmetry in \4)). 
The Hubbard- Strotonivich transformation used in our 
calculations is the spin-decomposition of Hirsch^i, which 
results in real Ising-like auxiliary fields. The phase prob- 
lem comes only from one-body hopping terms. Wc use a 
simple version of the phaseless approximatioi^i^ to con- 
strain 10) to a unique phase. At each step of propagation, 
the paths of the RWs are required to satisfy: 



>0, 



(4) 



where and \(f>') are the current and proposed po- 
sitions. The left-hand side is used in the importance 
samplingi^ii^i2^. Wc use the free-electron wave function 
as I^'t). Since this is an eigenfunction of the complex ki- 
netic energy terms of H, all the phase effect is absorbed 
in the deterministic one-body part. The condition on 
the RWs is equivalent to the original constrained path 
approximationi^, to which Eq. ([4]) reduces if = 0. The 
phase constraint in Eq. Q is the only approximation in 
our method. 

Since the approximation involves only the overall 
sign/phase of the many-body wave function, it is reason- 
able to expect that the results will be relatively insensi- 
tive to I^t)- Extensive benchmarks have shown this to 
be the case. The general approach has, in a variety of 
systemai^i^^i^i^^, given results among the most accurate 
that can be achieved presently from QMC. 

As a quantative measure in the current case, we com- 
pare e{n) in 3 X 3 Hubbard lattices (C/ = 4 and 8) 
between our method (QMC) and exact diagonalization 
(ED). At each density (both = Ni and the polar- 
ized case N-^ — Ni = 1, with A^j^ = 1,2,3,4), we calcu- 
late the ground-state energies for 1000 random values 
(identical in QMC and ED), average the results, and es- 
timate a statistical error bar. In the QMC results, the 
error bar is the combined statistical errors from the ran- 
dom distribution and the QMC sampling, although 
the latter is much smaller compared to the former in this 
system. The results are shown in Fig. [TJ The agreement 
between QMC and exact results is excellent. The relative 
error [cqmcC'T-) — eED('T-)]/|eED('T-)|, shown in the bottom 
panel, is essentially zero for {/ = 4 and is less than 1.5% 
for U = 8, across the entire density range. 



III. RESULTS 

A. Equation of state 

Our main energy results arc summarized in Fig.'s [2] 
and [3] In Fig. [21 the equation of state is presented for 
several lattice sizes and interaction strengths. For densi- 
ties n < 0.9, convergence of the averaged energy is rapid 
with respect to the set of random twists, and typically 
20 0's is sufficient. For densities closer to half- filling. 
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FIG. 2: (Color online) Ground state energy per site of the 2- 
D Hubbard model vs. density, for several interaction strengths 
and lattice sizes. Error bars are combined QMC and 0- 
integration statistical errors. As a result of TABC, curves are 
smooth and different lattice sizes are indistinguishable. The 
inset shows convergence to the thermodynamic limit with a 
magnified view. (To reduce clutter, only every fifth density is 
shown for each size.) It also illustrates the accuracy of the fit 
eflt(n) across the density range for the phase below Uc- 
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FIG. 3: (Color online) The hole energy ehih) vs. hole density 
h for interacting systems, derived from Fig. [2l A clear min- 
imum is seen for [/ > 4, at finite hole density /ic. The inset 
shows eh{h) for non- interacting Hubbard model calculated for 
lattices up to 40 x 40 with PBC. Note the kinks and the fiat 
part of the curves near half-filling. The magenta curve is for 
a 12 X 12 lattice (and the dashed line, 40 x 40) using TABC, 
which effectively eliminates the finite-size and shell effects. 



the energy has stronger fluctuations with 0. Further, 
the requirement on statistical accuracy is higher in this 
region, because the error bar on eh(h) is magnified by 
\/h (see Eq. In this case, the number of boundary 

conditions is increased (to 60-300). In each region, the 
same set of random values are used to help correlate 
the results at different densities. The main graph shows 
results from a Trotter time step At = 0.05; the fit below 
[Eq. ([5])] and results in the inset have been extrapolated 
to At = 0. Convergence to the thermodynamic limit is 
seen with all three lattice sizes in the main graph. As the 
inset shows, 12 x 12 and 16 x 16 are indistinguishable to 
within statistical errors (~ 10-3). 

In Fig. [31 the hole energy eh{h) derived from e(n) is 
plotted. The inset illustrates the large finite-size and 
shell effects under the usual PBC. Because of degenera- 
cies at the Fermi surface, the hole energy has kinks and 
is a constant below a finite hole concentration^. As the 
system size is increased, the Chih) curves show conver- 
gence, but only slowly. Indeed a false signal for PS is seen 
in the non-interacting systems. These features are re- 
moved by TABC, with which a smooth monotonie curve 
is obtained. Excellent convergence toward the thermo- 
dynamic limit is achieved with a 12 x 12 lattice. 

Interacting systems show similar behaviors: under 
PBC the same kinks appear in the e{n) vs. n curve o^^i'^° 
for the interaction strengths considered here. The combi- 
nation of CPMC and TABC leads to a dramatic improve- 
ment. The main panel of Fig. Oshows the hole energy for 
interacting systems. A clear minimum in eh{h) can be 
seen at a finite hole density in all cases when [/ > 4. At 
U = A. hcis ^ 0.07-0.1. As U is increased, the position of 
the minimum is seen to shift to the right, i.e., to a larger 



he- As U decreases to J7 = 2, eh{h) appears to decrease 
monotonically down to h ^ 0.014, the lowest doping in 
these lattices, although it cannot be completely ruled out 
a shallow (< 0.03 from e/i(0)) minimum exists within the 
statistical error bars. 

The energy results indicate that, near half-filling, the 
system phase-separates into a hole-free phase of density 
n = 1 and a phase at ric ^ 1 — he- Within a single phase, 
our results are expected to be at or near the thermody- 
namic limit. If the system is in a mixed state with two 
or more phases present, however, there are likely finite- 
size and/or interface effects. This appears to be the case 
from the data, where we see a minimum in the hole energy 
curves (as opposed to a flat region), as well as size varia- 
tions in eh{h) in the hole density range < /i < he- Sim- 
ilarly, if the system is in a spatially inhomogeneous spin 
or charge density wave state with very long wavelength 
modulations, for example a stripelike state with only one 
stripe in a lattice of linear dimension up to L 16, flnite- 
size effects would likely make it indistinguishable from a 
phase-separated state in our calculations. 

As a simple way to characterize the equation of state in 
the thermodynamic limit at n < ric, we fit the calculated 
e(n) on n S (0, 0.9) (size L > 12 only) to a 4-th order 
polynomial. For U — A this gives 

efit(n) = -4.004 71 + 3. 769 0.700 71^ + 0.091 n''. (5) 

Statistical errors in the fitted coefficients are 10"'^ to 
IQ-^. The inset in Fig. [5] shows the quality of the fit. 
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FIG. 4: (Color online) Spin-spin correlation function C(r) for 
a 12 X 12 Hubbard lattice with U = 4. Within the PS region, 
the system exhibits long-range AF correlation. The strength 
of the long-range correlation decreases with doping, and van- 
ishes at smaller densities. The inset shows a 16 x 16 lattice 
at i7 = 4, at a few selected densities near Uc'. n = 0.9688 
(blue square), 0.9453 (green diamond), 0.9297 (red empty di- 
amond), and 0.9063 (black empty square). To aid the eye, the 
absolute value |C(r)| is shown, along two separate directions. 
The behavior of the curves indicates the finite sizes of the AF 
phase in the periodic lattice. 



B. Spin-Spin Correlation 

At n = 1, the ground state is known to exhibit long- 
range antiferromagnetic (AF) orderi ^^i^^ Doping intro- 
duces frustration and tends to destroy the AF order. 
To see how this occurs and the relation to PS, we use 
the back-propagation tcchniqu o^^i^^ to calculate the spin- 
spin correlation function: 

^^^^ " ^ 5Z'^'^"j+"-.T - '^j+r,i)("j,T - '^j.i))' (6) 



where r is a vector on the lattice and (..) denotes expec- 
tation with respect to the ground state. The results for a 
12 X 12 lattice at [/ = 4, after twist-averaging, are shown 
in Fig. [31 AF order is evident at n = 1, as expected. 
Note that the magnitude of the long-range part is ~ 0.2, 
and double occupancy of t and j-electrons is significant, 
as the strength of the interaction U is moderate. The 
long-range order decays rapidly with n and, in the ho- 
mogeneous phase, only short range correlation remains. 
(The minimum of e^(/i) is around n = 0.9167 in 12 x 12.) 

A more quantitative picture can be seen from the spin 
structure factor: 5(q) = J2r C'(r)e"^ ''. When the system 
has AF order, S{q} will peak at {tt,tt). The calculated 
results are plotted in Fig. [5l as a function of n for three 
different lattice sizes. There is a marked difference be- 
tween the small and larger doping regions. Below a crit- 
ical density [n < Uc), S{'K,Tr) remains finite but is small 
and independent of lattice size, indicating the presence 
of short-range spin correlation but no long-range mag- 
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FIG. 5: (Color online) Spin structure factor at q = (■7i",7r) for 
three system sizes calculated ai U — 4. The lines are guides 
to the eye. The inset shows S{tv, tt) vs. lattice size at several 
densities (obtained by linear interpolation if the exact n is not 
available in the particular lattice). 



netic order. Beyond ric, 5'(7r, tt) increases quickly as n 
approaches 1. As the inset illustrates, at each density 
5(7r, tt) grows proportionally with system size, suggest- 
ing the presence of long-range AF order. 

We now further examine the spatial dependence of the 
spin correlation. From the Maxwell construction, the size 
of the AF region in a phase-separated system (n > ric) 
is A'af = (1 ~ h/hc) N. In our calculations, C(r) is av- 
eraged over imaginary-time and MC configurations. An 
AF cluster of linear dimension Zaf > L/2 should, due 
to "winding" around the periodic lattice, have a finite, 
constant tail |C(r)| beyond |r| ^ L — Iaf, while a smaller 
cluster should have a tail at zero beyond |r| ~ Zap- Our 
C(r) results are consistent with this. In 12 x 12, finite 
resolution gives only a handful of densities on the inter- 
val (nc,n), so Iaf is close to either L or 0, and we see 
long plateaus. The inset in Fig. |4] shows 16 x 16 lattices, 
focusing on several densities near ric- At n = 0.9688 
and 0.9453, Iaf > but the former (large Iaf) has 

a long flat tail while the latter shows a decline with |r| 
in the middle, indicating reduced contributions in the 
sum in Eq. ([6|). Similar effects are seen in the other pair 
{Iaf < L/2), with n ~ 0.9297 showing an extended in- 
termediate region in which |C(r)| is finite but decreasing, 
before the vanishing tail. 



IV. DISCUSSION AND CONCLUSION 

As we have discussed in Sec. II. B, our calculations use 
a non-perturbative, many-body QMC method. We re- 
turn again to the only approximation in the method, 
namely the phase constraint, to help further gauge its 
impact. Although the possibility of a systematic bias 
cannot be ruled out, every indication has been that our 
results are very accurate — including the quality of the 
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present data, the consistency between the energy and 
spin correlation results, and the extensive benchmarks 
to date. As mentioned, the constrained path approxi- 
mation has been tested (Refs. [l5l|2.j|3d and others) in 
various Hubbard systems under periodic or open bound- 
ary conditions. Accurate energy results are obtained. 
In realistic electronic systems, an approximation which 
is based on the same framework but which has to deal 
with a real two-body phase problem (as opposed to the 
non-stochastic one-body hopping phase here) has been 
benchmarked in molecules (Refs. Il6ll25ll26ll27l and others) 
against density-matrix renormalization group and quan- 
tum chemistry methods. Again the accuracy in the calcu- 
lated ground state enerfy is consistent with that of Fig.[T] 
In addition, several other factors in the present work 
provide more self-consistency checks and show the ro- 
bustness of the results. At n = 1 and L/ = 4, an exact en- 
ergy can be obtained with PBC: e(l) = -0.8618(2)i2i32, 
which is below our result: —0.8559(4)2^. Since our largest 
systematic error is expected to occur here (maximum n) , 
this suggests that the tendency for PS would, if anything, 
be underestimated by our energies. Under TABC the en- 
tire density range (including half-filling) is treated with 
the same approach. All calculations use the correspond- 
ing free-electron wave function as I^Pt)- An identical pro- 
cedure is applied which has no tuning or adjustable pa- 
rameters. Clearly the constraining |^'t) has no minimum 
in eh, but an unambiguous minimum emerges from the 



calculations. Neither does \^t) contain spin order, but 
the AF ordering appears and vanishes, consistently with 
the behavior of the energy. 

In summary, recent advances in QMC techniques have 
enabled us to determine the equation of state numerically 
in the 2-D Hubbard model at intermediate interactions. 
Our results show that, upon doping, the ground state 
separates into one phase with AF order (hole-free) and 
the rest without [iic ~ 0.92 for U ~ 4). (The nature of 
the spatially inhomogcneous state will require further in- 
vestigation, for example, the distinction between a phase- 
separated state in finite lattices and density waves with 
long wavelengths, as discussed in Sec. HI. A. More calcu- 
lations are on-going, which we plan to report in a future 
publication.) The size of the AF spin-density wave region 
vanishes at tIc, causing the system to lose long-range AF 
order. 
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